#BiocManager::install("clusterProfiler")
#install.packages("enrichplot")
#installed.packages("patchwork")
#BiocManager::install("SingleR")
#BiocManager::install("celldex")
#BiocManager::install("ensembldb")
#BiocManager::install('limma')
#BiocManager::install("DESeq2")
#organism = "org.Mm.eg.db"
#BiocManager::install(organism)
#BiocManager::install("pathview")
#install.packages("Seurat")
#install.packages("remotes")
#remotes::install_github("satijalab/seurat")
#BiocManager::install("rhdf5")
#install.packages("hdf5r")
#BiocManager::install("NeuCA")
#install.packages('devtools')
#install_github("ggjlab/scMCA")
#install.packages("pheatmap")
#BiocManager::install("scater")
#devtools::install_github("sqjin/CellChat")
#BiocManager::install("ComplexHeatmap")
In this step, we are loading the libraries which are required for analysis and plotting of the data.
library(Seurat)
library(devtools)
library(ggplot2)
library(patchwork)
library(dplyr)
library(clusterProfiler)
library(enrichplot)
# we use ggplot2 to add x axis labels (ex: ridgeplot)
library(stringr)
library(celldex)
library(ensembldb)
library(SingleR)
library(limma)
library(DESeq2)
library(org.Mm.eg.db)
library(rhdf5)
library(hdf5r)
library(pheatmap)
library(ComplexHeatmap)
library(CellChat)
library(patchwork)
options(stringsAsFactors = FALSE)
library(NMF)
library(ggalluvial)
my_cols <- c('B cells'='#F8766D','DC'='#39568CFF','Endothelial cells'='#CD9600','Epithelial cells'='#00C19A','Fibroblasts'= "#C77CFF",
'Macrophages'='#00A9FF','Monocytes'='#0CB702','Stromal cells'='#FF61CC', "ILC" = "grey", "Neutrophils" = "Darkred")
# Mouse 1
bcg_mouse_1 <- "/Volumes/Seagate_4TB/spatial_ntm_samples/demultiplexed/bcg_sample_round_one/outs/"
list.files(bcg_mouse_1)
## [1] "analysis" "cloupe.cloupe"
## [3] "filtered_feature_bc_matrix" "filtered_feature_bc_matrix.h5"
## [5] "metrics_summary.csv" "molecule_info.h5"
## [7] "possorted_genome_bam.bam" "possorted_genome_bam.bam.bai"
## [9] "probe_set.csv" "raw_feature_bc_matrix"
## [11] "raw_feature_bc_matrix.h5" "spatial"
## [13] "spatial_enrichment.csv" "web_summary.html"
bcg_mouse_1_data <- Load10X_Spatial(
bcg_mouse_1,
filename = "filtered_feature_bc_matrix.h5",
assay = "spatial",
slice = "bcg",
filter.matrix = TRUE,
to.upper = FALSE,
image = NULL)
# Mouse 2
bcg_mouse_2 <- "/Volumes/Seagate_4TB/spatial_ntm_samples/demultiplexed/bcg_sample_round_two/outs/"
list.files(bcg_mouse_2)
## [1] "analysis" "bcg_only_web_summary.html"
## [3] "cloupe.cloupe" "filtered_feature_bc_matrix"
## [5] "filtered_feature_bc_matrix.h5" "metrics_summary.csv"
## [7] "molecule_info.h5" "possorted_genome_bam.bam"
## [9] "possorted_genome_bam.bam.bai" "probe_set.csv"
## [11] "raw_feature_bc_matrix" "raw_feature_bc_matrix.h5"
## [13] "spatial" "spatial_enrichment.csv"
bcg_mouse_2_data <- Load10X_Spatial(
bcg_mouse_2,
filename = "filtered_feature_bc_matrix.h5",
assay = "spatial",
slice = "bcg",
filter.matrix = TRUE,
to.upper = FALSE,
image = NULL)
# Mouse 1
bcg_mouse_1_data <- SCTransform(bcg_mouse_1_data, assay = "spatial", verbose = FALSE)
# Mouse 2
bcg_mouse_2_data <- SCTransform(bcg_mouse_2_data, assay = "spatial", verbose = FALSE)
p1 <- SpatialFeaturePlot(bcg_mouse_1_data, features = "Cd19", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p2 <- SpatialFeaturePlot(bcg_mouse_1_data, features = "Igha", alpha = c(0.1, 1), pt.size.factor = 1.5)
p3 <- SpatialFeaturePlot(bcg_mouse_1_data, features = "Bcl6", alpha = c(0.1, 1), pt.size.factor = 1.5)
p4 <- SpatialFeaturePlot(bcg_mouse_1_data, features = "Lta", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p5 <- SpatialFeaturePlot(bcg_mouse_1_data, features = "Ltb", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p6 <- SpatialFeaturePlot(bcg_mouse_1_data, features = "Cd4", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p7 <- SpatialFeaturePlot(bcg_mouse_1_data, features = "Mki67", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p8 <- SpatialFeaturePlot(bcg_mouse_1_data, features = "Cxcr5", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8
p1 <- SpatialFeaturePlot(bcg_mouse_2_data, features = "Cd19", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p2 <- SpatialFeaturePlot(bcg_mouse_2_data, features = "Igha", alpha = c(0.1, 1), pt.size.factor = 1.5)
p3 <- SpatialFeaturePlot(bcg_mouse_2_data, features = "Bcl6", alpha = c(0.1, 1), pt.size.factor = 1.5)
p4 <- SpatialFeaturePlot(bcg_mouse_2_data, features = "Lta", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p5 <- SpatialFeaturePlot(bcg_mouse_2_data, features = "Ltb", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p6 <- SpatialFeaturePlot(bcg_mouse_2_data, features = "Cd4", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p7 <- SpatialFeaturePlot(bcg_mouse_2_data, features = "Mki67", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p8 <- SpatialFeaturePlot(bcg_mouse_2_data, features = "Cxcr5", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8
# Mouse 1
bcg_mouse_1_analyzed <- RunPCA(bcg_mouse_1_data, assay = "SCT", verbose = FALSE)
bcg_mouse_1_analyzed <- FindNeighbors(bcg_mouse_1_analyzed, reduction = "pca", dims = 1:5)
bcg_mouse_1_analyzed <- FindClusters(bcg_mouse_1_analyzed, verbose = FALSE, resolution = 0.8)
bcg_mouse_1_analyzed <- RunUMAP(bcg_mouse_1_analyzed, reduction = "pca", dims = 1:5)
# Mouse 2
bcg_mouse_2_analyzed <- RunPCA(bcg_mouse_2_data, assay = "SCT", verbose = FALSE)
bcg_mouse_2_analyzed <- FindNeighbors(bcg_mouse_2_analyzed, reduction = "pca", dims = 1:5)
bcg_mouse_2_analyzed <- FindClusters(bcg_mouse_2_analyzed, verbose = FALSE, resolution = 0.8)
bcg_mouse_2_analyzed <- RunUMAP(bcg_mouse_2_analyzed, reduction = "pca", dims = 1:5)
# First convert the seurat object to single cell experiment object (otherwise the SingleR pipeline will not work)
# Mouse 1
bcg_mouse_1_analyzed.sce <- as.SingleCellExperiment(bcg_mouse_1_analyzed)
# Mouse 2
bcg_mouse_2_analyzed.sce <- as.SingleCellExperiment(bcg_mouse_2_analyzed)
#create reference data
ref.data <- celldex::ImmGenData()
# run singleR pipeline to find the cell types
## Mouse 1
cell_types_mouse_1 <- SingleR(test=bcg_mouse_1_analyzed.sce, assay.type.test=1,
ref=ref.data, labels=ref.data$label.main)
## Mouse 2
cell_types_mouse_2 <- SingleR(test=bcg_mouse_2_analyzed.sce, assay.type.test=1,
ref=ref.data, labels=ref.data$label.main)
# view cell types
# Mouse 1
table(cell_types_mouse_1$labels)
##
## B cells DC Endothelial cells Epithelial cells
## 48 43 1244 264
## Fibroblasts Macrophages Monocytes Neutrophils
## 437 1052 20 2
## Stromal cells
## 208
#checking cell types Mouse 1
plotScoreHeatmap(cell_types_mouse_1)
plotDeltaDistribution(cell_types_mouse_1)
summary(is.na(cell_types_mouse_1$pruned.labels))
## Mode FALSE
## logical 3318
# Mouse 2
table(cell_types_mouse_2$labels)
##
## B cells DC Endothelial cells Epithelial cells
## 130 62 1027 175
## Fibroblasts Macrophages Monocytes Neutrophils
## 335 707 5 1
## Stromal cells
## 258
#checking cell types Mouse 2
plotScoreHeatmap(cell_types_mouse_2)
plotDeltaDistribution(cell_types_mouse_2)
summary(is.na(cell_types_mouse_2$pruned.labels))
## Mode FALSE TRUE
## logical 2694 6
# Mouse 1
bcg_mouse_1_analyzed[["ref.data"]] <- cell_types_mouse_1$labels
# Mouse 2
bcg_mouse_2_analyzed[["ref.data"]] <- cell_types_mouse_2$labels
p1 <- DimPlot(bcg_mouse_1_analyzed, group.by = c("ref.data"), reduction = "umap",label = F, pt.size = 1.5, label.size = 5, cols = my_cols)
p2 <- SpatialDimPlot(bcg_mouse_1_analyzed, group.by = c("ref.data"), label = FALSE, label.size = 5, cols = my_cols)
p1 + p2
p1 <- DimPlot(bcg_mouse_2_analyzed, group.by = c("ref.data"), reduction = "umap",label = F, pt.size = 1.5, label.size = 5, cols = my_cols)
p2 <- SpatialDimPlot(bcg_mouse_2_analyzed, group.by = c("ref.data"), label = FALSE, label.size = 5, cols = my_cols)
p1 + p2
VlnPlot(bcg_mouse_1_analyzed, group.by = c("ref.data"), features = c("Cxcr5", "Mki67", "Ltb", "Cd38", "Ccr6", "Tnfrsf13c"), cols = my_cols)
VlnPlot(bcg_mouse_2_analyzed, group.by = c("ref.data"), features = c("Cxcr5", "Mki67", "Ltb", "Cd38", "Ccr6", "Tnfrsf13c"), cols = my_cols)
B_cell_mouse_1 <- FindMarkers(bcg_mouse_1_analyzed, group.by = "ref.data", ident.1 = "B cells", min.pct = 0.25)
head(B_cell_mouse_1, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Cr2 8.320332e-60 0.9384812 0.604 0.052 1.291315e-55
## Ms4a1 1.615581e-56 0.8401004 0.604 0.055 2.507382e-52
## Ighd 4.694334e-53 1.3861847 0.812 0.120 7.285606e-49
## Bank1 8.505742e-52 0.7844861 0.646 0.068 1.320091e-47
## Cd19 4.140845e-51 1.4892818 0.917 0.167 6.426591e-47
## Il9r 3.138481e-50 0.4293062 0.354 0.020 4.870923e-46
## Blnk 3.244900e-45 1.0068607 0.812 0.132 5.036085e-41
## Tnfrsf13c 1.132859e-43 1.1269492 0.750 0.119 1.758198e-39
## Spib 2.486053e-39 1.2094113 0.708 0.117 3.858354e-35
## Pla2g2d 2.630705e-38 0.7144449 0.479 0.050 4.082853e-34
write.csv(B_cell_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/B_cell_mouse_1.csv")
DC_mouse_1 <- FindMarkers(bcg_mouse_1_analyzed, group.by = "ref.data", ident.1 = "DC", min.pct = 0.25)
head(DC_mouse_1, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Fam241a 1.741273e-09 0.2699921 0.279 0.058 2.702456e-05
## Slfn9 6.331009e-09 0.2754009 0.256 0.053 9.825725e-05
## Fmo2 9.853839e-09 -0.9743768 0.349 0.690 1.529316e-04
## Epas1 1.523917e-08 -0.8608185 0.605 0.896 2.365120e-04
## Plxnc1 3.464226e-08 0.4707331 0.535 0.199 5.376479e-04
## Igfbp2 7.898749e-08 -0.8938734 0.558 0.867 1.225886e-03
## Slpi 4.991795e-07 0.4812667 0.953 0.868 7.747265e-03
## Aqp1 9.984770e-07 -0.7826018 0.233 0.595 1.549636e-02
## Ltb 1.031783e-06 0.3426673 0.465 0.175 1.601327e-02
## Sftpc 1.082770e-06 -0.2978901 1.000 0.999 1.680459e-02
write.csv(DC_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/DC_mouse_1.csv")
Stromal_cells_mouse_1 <- FindMarkers(bcg_mouse_1_analyzed, group.by = "ref.data", ident.1 = "Stromal cells", min.pct = 0.25)
head(Stromal_cells_mouse_1, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Igkc 2.549004e-17 -0.5568178 1.000 0.999 3.956055e-13
## Igha 1.333650e-16 -0.5432519 1.000 0.997 2.069825e-12
## Gsn 1.929346e-15 0.4734754 0.933 0.764 2.994344e-11
## Jchain 5.250575e-15 -0.6192559 0.654 0.849 8.148893e-11
## Tnfaip2 2.023024e-13 -0.4416083 0.933 0.964 3.139733e-09
## Epas1 3.849126e-12 0.3735630 0.995 0.885 5.973844e-08
## Inmt 4.073284e-12 0.4501445 0.981 0.841 6.321736e-08
## Ighj4 1.951249e-11 -0.5177593 0.630 0.801 3.028338e-07
## Timp3 2.877549e-11 0.3253112 0.995 0.967 4.465955e-07
## Cldn5 2.079806e-10 0.3490426 0.952 0.835 3.227860e-06
write.csv(Stromal_cells_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Stromal_cells_mouse_1.csv")
Endothelial_cells_mouse_1 <- FindMarkers(bcg_mouse_1_analyzed, group.by = "ref.data", ident.1 = "Endothelial cells", min.pct = 0.25)
head(Endothelial_cells_mouse_1, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Igfbp2 2.960434e-64 0.6943825 0.940 0.817 4.594593e-60
## Epas1 4.555557e-63 0.5417082 0.957 0.853 7.070224e-59
## Cldn5 7.254164e-62 0.5775664 0.934 0.788 1.125846e-57
## Tns1 3.981420e-47 0.4378691 0.965 0.883 6.179164e-43
## Ace 1.107643e-46 0.4589540 0.747 0.528 1.719061e-42
## Spock2 2.544447e-44 0.4939606 0.722 0.540 3.948982e-40
## Pakap.1 6.648450e-44 0.4067602 0.920 0.843 1.031840e-39
## Rasip1 6.994426e-43 0.4084629 0.716 0.498 1.085535e-38
## Timp3 9.312203e-43 0.3806280 0.994 0.953 1.445254e-38
## Clec14a 4.269361e-42 0.3756253 0.589 0.358 6.626049e-38
write.csv(Endothelial_cells_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Endothelial_cells_mouse_1.csv")
Epithelial_cell_mouse_1 <- FindMarkers(bcg_mouse_1_analyzed, group.by = "ref.data", ident.1 = "Epithelial cells", min.pct = 0.25)
head(Epithelial_cell_mouse_1, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Inmt 3.735506e-25 -0.8116861 0.678 0.865 5.797505e-21
## Tns1 1.355276e-21 -0.5743788 0.818 0.922 2.103388e-17
## Igha 2.119608e-20 0.3895227 1.000 0.997 3.289632e-16
## Igfbp2 2.733816e-19 -0.7586417 0.712 0.876 4.242883e-15
## Cldn5 2.164530e-18 -0.6047003 0.708 0.854 3.359351e-14
## Epas1 4.393914e-18 -0.5764613 0.777 0.902 6.819354e-14
## Timp3 1.315169e-17 -0.4704262 0.928 0.972 2.041142e-13
## Pakap.1 3.745781e-17 -0.4711075 0.777 0.880 5.813452e-13
## Igkc 4.927450e-17 0.3794999 1.000 0.999 7.647402e-13
## Cd93 7.279343e-17 -0.4278941 0.178 0.446 1.129754e-12
write.csv(Epithelial_cell_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Epithelial_cell_mouse_1.csv")
Fibroblasts_mouse_1 <- FindMarkers(bcg_mouse_1_analyzed, group.by = "ref.data", ident.1 = "Fibroblasts", min.pct = 0.25)
head(Fibroblasts_mouse_1, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Inmt 5.509154e-26 0.4909881 0.966 0.832 8.550208e-22
## Gsn 1.321252e-25 0.4882381 0.906 0.754 2.050583e-21
## Ltbp4 1.503307e-17 0.3913844 0.654 0.474 2.333132e-13
## Clec3b 1.204074e-16 0.4400577 0.613 0.437 1.868723e-12
## Nbl1 2.244220e-15 0.3431696 0.883 0.774 3.483029e-11
## Fmo2 1.080310e-14 0.3284504 0.808 0.667 1.676641e-10
## Rarres2 2.062960e-14 0.4293750 0.751 0.625 3.201714e-10
## Timp3 1.157632e-13 0.2623184 0.993 0.965 1.796645e-09
## Gbp2b 2.138363e-13 -0.2928538 0.924 0.953 3.318739e-09
## Ogn 5.905678e-13 0.2581432 0.384 0.234 9.165613e-09
write.csv(Fibroblasts_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Fibroblasts_mouse_1.csv")
Macrophages_mouse_1 <- FindMarkers(bcg_mouse_1_analyzed, group.by = "ref.data", ident.1 = "Macrophages", min.pct = 0.25)
head(Macrophages_mouse_1, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Igfbp2 2.890373e-66 -0.7985205 0.787 0.898 4.485859e-62
## Epas1 1.058095e-59 -0.5990300 0.828 0.922 1.642164e-55
## Igha 7.269746e-58 0.4226417 0.997 0.997 1.128265e-53
## Inmt 1.475707e-55 -0.7200217 0.775 0.885 2.290297e-51
## Cldn5 5.228715e-52 -0.6074341 0.761 0.880 8.114965e-48
## Igkc 3.122946e-49 0.3467899 0.998 0.999 4.846812e-45
## Ctss 3.077046e-48 0.4188831 0.998 0.981 4.775575e-44
## Timp3 8.672736e-48 -0.4636797 0.952 0.976 1.346009e-43
## Tnfaip2 1.139208e-45 0.3668298 0.993 0.948 1.768051e-41
## Tns1 3.488936e-43 -0.4879362 0.872 0.934 5.414828e-39
write.csv(Macrophages_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Macrophages_mouse_1.csv")
Monocytes_mouse_1 <- FindMarkers(bcg_mouse_1_analyzed, group.by = "ref.data", ident.1 = "Macrophages", min.pct = 0.25)
head(Monocytes_mouse_1, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Igfbp2 2.890373e-66 -0.7985205 0.787 0.898 4.485859e-62
## Epas1 1.058095e-59 -0.5990300 0.828 0.922 1.642164e-55
## Igha 7.269746e-58 0.4226417 0.997 0.997 1.128265e-53
## Inmt 1.475707e-55 -0.7200217 0.775 0.885 2.290297e-51
## Cldn5 5.228715e-52 -0.6074341 0.761 0.880 8.114965e-48
## Igkc 3.122946e-49 0.3467899 0.998 0.999 4.846812e-45
## Ctss 3.077046e-48 0.4188831 0.998 0.981 4.775575e-44
## Timp3 8.672736e-48 -0.4636797 0.952 0.976 1.346009e-43
## Tnfaip2 1.139208e-45 0.3668298 0.993 0.948 1.768051e-41
## Tns1 3.488936e-43 -0.4879362 0.872 0.934 5.414828e-39
write.csv(Monocytes_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Monocytes_mouse_1.csv")
B_cell_mouse_2 <- FindMarkers(bcg_mouse_2_analyzed, group.by = "ref.data", ident.1 = "B cells", min.pct = 0.25)
head(B_cell_mouse_2, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Cxcr5 1.795952e-103 0.8833483 0.638 0.072 2.455426e-99
## Tnfrsf13c 1.976759e-98 1.6255464 0.877 0.182 2.702624e-94
## Cr2 2.568990e-98 1.1656010 0.623 0.075 3.512323e-94
## Cd19 2.894342e-97 1.6119812 0.885 0.195 3.957144e-93
## Iglc3 2.116963e-92 1.2005717 0.862 0.178 2.894312e-88
## Ighd 3.359851e-88 1.4803253 0.808 0.165 4.593588e-84
## Fcrla 1.286344e-85 0.9721240 0.677 0.105 1.758689e-81
## Bank1 1.955197e-84 1.0070053 0.692 0.112 2.673146e-80
## Pax5 2.727178e-83 0.9189419 0.662 0.102 3.728597e-79
## Ly6d 1.258899e-81 0.5611298 0.515 0.056 1.721167e-77
write.csv(B_cell_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/B_cell_mouse_2.csv")
DC_mouse_2 <- FindMarkers(bcg_mouse_2_analyzed, group.by = "ref.data", ident.1 = "DC", min.pct = 0.25)
head(DC_mouse_2, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Itk 1.126374e-24 0.6518767 0.806 0.240 1.539978e-20
## Sftpc 2.294687e-22 -1.2702568 1.000 1.000 3.137296e-18
## Pak1 2.942115e-20 0.5473896 0.726 0.218 4.022460e-16
## Clic5 4.911070e-20 -1.2458544 0.694 0.956 6.714415e-16
## C3 5.223048e-20 0.8231177 1.000 0.997 7.140951e-16
## Cd3d 5.577503e-20 0.8201085 0.839 0.343 7.625563e-16
## Ccr7 1.958530e-19 0.6471925 0.629 0.187 2.677702e-15
## Ager 5.354944e-19 -1.1661539 0.903 0.973 7.321280e-15
## H2-Eb1 7.763557e-19 0.6541973 1.000 0.999 1.061433e-14
## Vegfa 1.098108e-18 -1.1992835 0.629 0.932 1.501333e-14
write.csv(DC_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/DC_mouse_2.csv")
Stromal_cells_mouse_2 <- FindMarkers(bcg_mouse_2_analyzed, group.by = "ref.data", ident.1 = "Stromal cells", min.pct = 0.25)
head(Stromal_cells_mouse_2, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Gsn 1.386259e-28 0.6018293 0.996 0.905 1.895294e-24
## Tnfaip2 4.346716e-23 -0.7090419 0.953 0.984 5.942831e-19
## Psap 1.665689e-21 -0.4323633 1.000 1.000 2.277331e-17
## Myl9 1.259329e-19 0.4828020 0.477 0.239 1.721755e-15
## H2-DMb2 1.679277e-19 -0.6117161 0.946 0.977 2.295907e-15
## Ctss 4.810568e-19 -0.5572010 0.992 0.995 6.577009e-15
## Ly6e 5.775611e-19 -0.4271559 0.984 0.998 7.896415e-15
## B2m 1.406355e-18 -0.3576475 0.996 0.998 1.922769e-14
## Cyp2f2 7.401773e-18 0.8529229 1.000 0.941 1.011970e-13
## Myh11 3.177795e-17 0.5537961 0.473 0.249 4.344681e-13
write.csv(Stromal_cells_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Stromal_cells_mouse_2.csv")
Endothelial_cells_mouse_2 <- FindMarkers(bcg_mouse_2_analyzed, group.by = "ref.data", ident.1 = "Endothelial cells", min.pct = 0.25)
head(Endothelial_cells_mouse_2, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Cldn5 8.891607e-143 0.9344905 0.985 0.738 1.215661e-138
## Igfbp2 4.848456e-140 1.0061615 0.983 0.708 6.628810e-136
## Tns1 1.456044e-130 0.7875292 1.000 0.882 1.990704e-126
## Epas1 2.648369e-124 0.7706319 0.997 0.854 3.620850e-120
## Inmt 2.626647e-108 0.8388433 0.979 0.682 3.591152e-104
## Timp3 2.201202e-102 0.6322175 0.998 0.929 3.009484e-98
## Clic5 9.961688e-101 0.5969119 0.999 0.921 1.361962e-96
## Pakap.1 1.276464e-99 0.6482757 0.990 0.892 1.745182e-95
## Spock2 1.998597e-97 0.7174474 0.937 0.628 2.732482e-93
## Acer2 2.689486e-93 0.6562824 0.951 0.763 3.677065e-89
write.csv(Endothelial_cells_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Endothelial_cells_mouse_2.csv")
Epithelial_cell_mouse_2 <- FindMarkers(bcg_mouse_2_analyzed, group.by = "ref.data", ident.1 = "Epithelial cells", min.pct = 0.25)
head(Epithelial_cell_mouse_2, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Igfbp2 4.276103e-42 -1.7397435 0.474 0.836 5.846287e-38
## Inmt 4.796397e-41 -1.6731033 0.474 0.817 6.557634e-37
## Spint2 3.086740e-40 0.7716912 1.000 0.972 4.220190e-36
## Ccdc96 9.863237e-40 0.4891312 0.326 0.059 1.348502e-35
## Wfdc2 3.877575e-38 1.3042103 0.943 0.806 5.301421e-34
## Tns1 2.295413e-37 -1.1410230 0.783 0.937 3.138288e-33
## Cdhr4 3.244769e-37 0.7131494 0.400 0.096 4.436249e-33
## Erich3 6.609781e-37 0.4512833 0.320 0.060 9.036892e-33
## 1700001C02Rik 1.067681e-36 0.3713457 0.303 0.053 1.459733e-32
## Cldn5 1.185186e-36 -1.3527028 0.583 0.849 1.620386e-32
write.csv(Epithelial_cell_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Epithelial_cell_mouse_2.csv")
Fibroblasts_mouse_2 <- FindMarkers(bcg_mouse_2_analyzed, group.by = "ref.data", ident.1 = "Fibroblasts", min.pct = 0.25)
head(Fibroblasts_mouse_2, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Cst3 4.019251e-33 0.4103062 1.000 1.000 5.495120e-29
## Gsn 6.675412e-26 0.5246076 0.976 0.905 9.126623e-22
## Clec3b 5.877289e-24 0.5597167 0.887 0.725 8.035430e-20
## Mgp 1.809036e-23 0.3799454 1.000 0.992 2.473314e-19
## Ltbp4 9.015487e-22 0.5156064 0.866 0.682 1.232597e-17
## Inmt 1.181286e-21 0.5020674 0.961 0.771 1.615054e-17
## Eln 6.146593e-17 0.4589563 0.815 0.668 8.403622e-13
## H2-DMb2 2.324752e-16 -0.4986663 0.967 0.975 3.178401e-12
## H2-Ab1 1.293291e-15 -0.3523212 0.988 0.991 1.768188e-11
## Ctss 1.964660e-15 -0.4255406 0.997 0.995 2.686083e-11
write.csv(Fibroblasts_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Fibroblasts_mouse_2.csv")
Macrophages_mouse_2 <- FindMarkers(bcg_mouse_2_analyzed, group.by = "ref.data", ident.1 = "Macrophages", min.pct = 0.25)
head(Macrophages_mouse_2, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Ctss 1.040821e-95 0.7567875 1.000 0.993 1.423011e-91
## Psap 3.925294e-89 0.5558869 1.000 1.000 5.366662e-85
## Epas1 2.949012e-88 -0.8781662 0.827 0.937 4.031889e-84
## Inmt 6.472508e-85 -1.1013762 0.611 0.860 8.849213e-81
## Gbp2b 9.581003e-81 0.6528410 0.993 0.962 1.309915e-76
## Tnfaip2 2.300836e-79 0.7575139 0.997 0.976 3.145703e-75
## Cldn5 4.546197e-76 -0.9452875 0.699 0.879 6.215561e-72
## Ctsc 1.261996e-71 0.6238895 0.997 0.957 1.725401e-67
## Nos2 1.006443e-70 0.9034125 0.755 0.441 1.376009e-66
## Igfbp2 4.797133e-70 -1.0009458 0.666 0.865 6.558640e-66
write.csv(Macrophages_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Macrophages_mouse_2.csv")
Monocytes_mouse_2 <- FindMarkers(bcg_mouse_2_analyzed, group.by = "ref.data", ident.1 = "Monocytes", min.pct = 0.25)
head(Monocytes_mouse_2, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Epha8 1.726519e-32 0.4779517 0.4 0.005 2.360496e-28
## Zcchc18 1.701735e-23 0.4752914 0.4 0.007 2.326611e-19
## Nxpe4 3.777130e-19 0.4726359 0.4 0.009 5.164092e-15
## Saa3 1.956962e-12 0.4641713 0.4 0.014 2.675559e-08
## Cacng8 3.698153e-11 0.4610098 0.4 0.016 5.056115e-07
## Tmem69 5.716732e-10 0.4589060 0.4 0.019 7.815916e-06
## Cd5 8.543622e-10 0.4583805 0.4 0.019 1.168084e-05
## Paip1 3.760272e-09 0.4557560 0.4 0.020 5.141044e-05
## Olfm2 4.327556e-09 0.4494765 0.4 0.020 5.916634e-05
## 1810058I24Rik 1.008572e-08 0.4536598 0.4 0.022 1.378920e-04
write.csv(Monocytes_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Monocytes_mouse_2.csv")
Macrophage_cell_deseq <- FindMarkers(bcg_mouse_1_analyzed, group.by = "ref.data", ident.1 = "Macrophages", test.use = "DESeq2")
write.csv(Macrophage_cell_deseq, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Macrophage_cell_deseq.csv")
bcg_cellchat <- createCellChat(object = bcg_mouse_2_analyzed, group.by = "ref.data")
## The `data` slot in the default assay is used. The default assay is SCT
## The `meta.data` slot in the Seurat object is used as cell meta information
## The cell groups used for CellChat analysis are B cells DC Endothelial cells Epithelial cells Fibroblasts Macrophages Monocytes Neutrophils Stromal cells
CellChatDB <- CellChatDB.mouse
showDatabaseCategory(CellChatDB)
CellChatDB.use <- CellChatDB
bcg_cellchat@DB <- CellChatDB.use
bcg_cellchat <- subsetData(bcg_cellchat)
bcg_cellchat <- identifyOverExpressedGenes(bcg_cellchat)
bcg_cellchat <- identifyOverExpressedInteractions(bcg_cellchat)
bcg_cellchat <- computeCommunProb(bcg_cellchat)
bcg_cellchat <- filterCommunication(bcg_cellchat, min.cells = 10)
## The cell-cell communication related with the following cell groups are excluded due to the few number of cells: Monocytes Neutrophils
bcg_cellchat <- computeCommunProbPathway(bcg_cellchat)
bcg_cellchat <- aggregateNet(bcg_cellchat)
groupSize <- as.numeric(table(bcg_cellchat@idents))
par(mfrow = c(1,2), xpd=TRUE)
netVisual_circle(bcg_cellchat@net$count, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Number of interactions")
netVisual_circle(bcg_cellchat@net$weight, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Interaction weights/strength")
# See interaction of each cell type with the another
mat1 <- bcg_cellchat@net$weight
par(mfrow = c(3,4), xpd=TRUE)
for (i in 1:nrow(mat1)) {
mat2 <- matrix(0, nrow = nrow(mat1), ncol = ncol(mat1), dimnames = dimnames(mat1))
mat2[i, ] <- mat1[i, ]
netVisual_circle(mat2, vertex.weight = groupSize, weight.scale = T, edge.weight.max = max(mat1), title.name = rownames(mat1)[i])
}
##Visualize signalling pathways in each cell
# Identify signaling pathways showing significant communications
bcg_cellchat@netP$pathways
## [1] "COLLAGEN" "APP" "GALECTIN" "MHC-II" "FN1"
## [6] "THBS" "MIF" "LAMININ" "GAS" "VEGF"
## [11] "COMPLEMENT" "SEMA4" "SEMA3" "JAM" "ICAM"
## [16] "CCL" "ITGAL-ITGB2" "SPP1" "CD22" "CD45"
## [21] "PDGF" "TGFb" "CXCL" "GRN" "CHEMERIN"
## [26] "PARs" "TENASCIN" "PECAM1" "WNT" "TNF"
## [31] "NOTCH" "CD52" "PD-L1" "ESAM" "SEMA7"
## [36] "CDH" "EPHA" "IL16" "THY1" "ncWNT"
## [41] "IGF" "APRIL" "ICOS" "TWEAK" "PERIOSTIN"
## [46] "CDH5" "PDL2"
# Evaluate one of the pathways
pathways.show <- c("MHC-II")
vertex.receiver = seq(2,5) # a numeric vector.
?netVisual_aggregate(bcg_cellchat, signaling = pathways.show, vertex.receiver = vertex.receiver)
# Chord diagram
par(mfrow=c(1,1))
netVisual_aggregate(bcg_cellchat, signaling = pathways.show, layout = "chord")
# Heatmap
par(mfrow=c(1,1))
netVisual_heatmap(bcg_cellchat, signaling = pathways.show, color.heatmap = "Reds")
netAnalysis_contribution(bcg_cellchat, signaling = pathways.show)
pairLR.MHC <- extractEnrichedLR(bcg_cellchat, signaling = pathways.show, geneLR.return = FALSE)
LR.show <- pairLR.MHC[1,] # show one ligand-receptor pair
# Hierarchy plot
vertex.receiver = seq(1,4) # a numeric vector
netVisual_individual(bcg_cellchat, signaling = pathways.show, pairLR.use = LR.show, vertex.receiver = vertex.receiver)
## [[1]]
# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
netVisual_bubble(bcg_cellchat, sources.use = 1, targets.use = c(2:8), remove.isolate = FALSE)
# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
# show all the interactions sending from Inflam.FIB
netVisual_chord_gene(bcg_cellchat, sources.use = 1, targets.use = c(2:8), lab.cex = 0.5,legend.pos.y = 30)
bcg_cellchat <- netAnalysis_computeCentrality(bcg_cellchat, slot.name = "netP")
# the slot 'netP' means the inferred intercellular communication network of signaling pathways
# Visualize the computed centrality scores using heatmap, allowing ready identification of major signaling roles of cell groups
netAnalysis_signalingRole_network(bcg_cellchat, signaling = pathways.show, width = 8, height = 2.5, font.size = 10)
ht1 <- netAnalysis_signalingRole_heatmap(bcg_cellchat, pattern = "outgoing")
ht2 <- netAnalysis_signalingRole_heatmap(bcg_cellchat, pattern = "incoming")
ht1 + ht2
selectK(bcg_cellchat, pattern = "outgoing")
nPatterns = 3
bcg_cellchat <- identifyCommunicationPatterns(bcg_cellchat, pattern = "outgoing", k = nPatterns)
# river plot
netAnalysis_river(bcg_cellchat, pattern = "outgoing")
# dot plot
netAnalysis_dot(bcg_cellchat, pattern = "outgoing")
selectK(bcg_cellchat, pattern = "incoming")
nPatterns = 3
bcg_cellchat <- identifyCommunicationPatterns(bcg_cellchat, pattern = "incoming", k = nPatterns)
# river plot
netAnalysis_river(bcg_cellchat, pattern = "incoming")
# dot plot
netAnalysis_dot(bcg_cellchat, pattern = "incoming")
organism = org.Mm.eg.db
# reading in data from deseq2
df = read.csv("/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Macrophage_cell_deseq.csv", header=TRUE)
# we want the log2 fold change
original_gene_list <- df$avg_log2FC
# name the vector
names(original_gene_list) <- df$X
# omit any NA values
gene_list<-na.omit(original_gene_list)
# sort the list in decreasing order (required for clusterProfiler)
gene_list = sort(gene_list, decreasing = TRUE)
#gene set enrichment function for B cell cluster
gse <- gseGO(geneList=gene_list,
ont ="ALL",
keyType = "SYMBOL",
minGSSize = 3,
maxGSSize = 800,
pvalueCutoff = 0.05,
verbose = TRUE,
OrgDb = organism,
pAdjustMethod = "none",
eps = 0.0)
dotplot(gse, showCategory=10, split=".sign") + facet_grid(.~.sign)